setwd("E:\\5hmc_file\\2_5hmc_yjp_bam\\ASM\\bayes_pvalue_beta0")
library(tidyr)
library(dplyr)
group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39")
sel1=paste(group1,".bayes_p.txt",sep="")
rt=read.csv("pvalue_sig_statis.csv",header = T)
id=as.character(rt[rt$twins>0,]$unitID) #选出至少1对twins都是偏移的位点
i=1
for(i in 2:length(sel1)){
file=read.table(sel1[i],head=T,sep="\t")
file$unitID=paste(file$chrom,file$position,sep=":")
file=file[file$normal_bayes_pvalue<0.05&file$tumor_bayes_pvalue<0.05,]
file1=select(file,unitID,normal_bayes_beta0,tumor_bayes_beta0)
file1$normal_group1=ifelse(file1$normal_bayes_beta0>0,"normal.Up","normal.Down")
file1$tumor_group1=ifelse(file1$tumor_bayes_beta0>0,"tumor.Up","tumor.Down")
file1$group1=paste(file1$normal_group1,file1$tumor_group1,sep="_")
test=data.frame(table(file1$group1))
testnew=data.frame(matrix(NA,2,2,byrow = T))
names(testnew)=c("Var1","Freq")
testnew[1,1]="same_direction"
testnew[1,2]=test[1,]$Freq+test[4,]$Freq
testnew[2,1]="oppo_direction"
testnew[2,2]=test[2,]$Freq+test[3,]$Freq
test=rbind(test,testnew)
names(test)=c("Var1",group1[i])
fn=paste("../20201103/",group1[i],".pvalue.both.sig.twins.txt",sep="")
write.table(file1,fn,quote=F,row.names = F,sep="\t")
rt1=merge(rt1,test,by="Var1")
}
write.csv(rt1,"../20201103/pvalue.both.sig.twins.statis.csv",quote = F,row.names = F)
